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I . INTRODUCTION 



For many analysis problems in turbomachinery, an assumption 
of inviscid flow provides sufficiently accurate results for 
design or development tasks. At the same time the assumption 
or irrotational flow will limit application to low speed shock- 
less cases. The Euler equations describe inviscid rotational 
flow and therefore allow the modeling of transonic or super- 
sonic flows with shock waves and rotational flow fields behind 
strong waves. Solutions of the Euler equations are therefore 
important for many applications. 

A computer program was developed which solves numerically 
two-dimensional Euler equations using the Godunov method 
(Ref. 1) or Higher Order Extension of the Godunov method 
(Ref. 2) . The program was intended primarily for turbomachinery 
applications, but with modification of the boundary conditions, 
application to any problem which is modeled by the Euler equa- 
tions is possible. 

The purpose of this report is to give instructions to 
potential users of the present computer code and to outline 
the basic assumptions made in the mathematical model and 
numerical solution. 
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II. GOVERNING EQUATIONS 



The unsteady two-dimensional Euler equations can 
be written in conservation law form as 



9a 95 3c _ n 

3t 9x 3y 



( 1 ) 
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and 



e = p(e + — — y — ) = energy of the unit of volume 



e = Y y - 1) p ~ i nterna l energy 
p = density 

U and V are the velocity components corresponding to 
the x and y coordinates 

p = pressure 

Y = specific heat ratio 

Equation (1) is non-dimensionalized using a refer- 
ence length p Q , a reference pressure x and a reference 
density p Q . These reference quantities are in princi- 
ple arbitrary and were selected as follows 

x- = relative length of the geometrical cord 

p ,= 101350.0 -EL 
m 

P Q = 1.2 kg/m^ 
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The reference quantities for time, t and velocity 
not independent and are determined from 



a are 
o 




1 ; 



( 3 ) 



then 



a 



o 



J 




t = x /a 
o o / o ; 



non-dimensionalize the variables in Eq. (1) as follows; 



When these quantities are inserted in Eq. (1) , the flow 
equations, in terms of non-dimensional variables, are 
identical in form to those in terms of dimensional 
variables. The subscript "n" may then be dropped without 
confusion and only non-dimensional quantities need be 
referred to in the remaining sections of this report. 

We will solve Eq. (1) in physical space, so no further 
transformation of the equation is required. 





x = x _ x _ 
o n 




( 4 ) 
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III. BOUNDARY CONDITIONS 



The computational domain illustrated in Figure 1 
is bounded by open surfaces (flow may cross these sur- 
faces) , solid surfaces and inflow-outflow boundaries. 
Formulation and implementation of boundary conditions at 
each surface are discussed in this section. 

3.1. SOLID WALL BOUNDARY CONDITION 



The solid wall boundary condition appropriate for 
inviscid flow computations is that of zero mass flux 
through the surface. This condition is difficult to im- 
plement uniquely for the Euler equations. 

In the program developed, the user has a choice of 



two types of boundary conditions on the solid walls as 
follows ; 



a) Solution of the Riemann problem between the 
point of the domain of integration and its 
mirror image in the direction normal to the 
wall with the same parameters but, with an 
inverse sign of velocity component normal to 
the wall. (See Fig. 2) . 



b) 



The wall pressure is derived from the momen- 
tum equation 




9n R 

s 



(5) 



where 



n - indicates the direction normal to the 
surface 
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Up = velocity parallel to the surface 

R = radius of curvature of the surface 
s 

If we add the flow tangency condition, we have a 
uniquely defined boundary condition at the wall. The 
solid wall boundary conditions are imposed at the cascade 
surface, which corresponds to the segments 3-5 and 4-6 
in the computational domain shown in Figure 1. 

3.2. PERIODIC BOUNDARIES 

The periodic flow boundary conditions are imposed 
on the segments 1-3, 5-7, 2-4, and 6-8 in Fig. 1. 

By virtue of the periodicity of the flow with respect to 
y, the flow parameters over segments 1-3 and 5-7 are the 
same as those over 2-4 and 6-8. Hence, the segments 1-3, 
5-7, 2-4, and 6-8 are essentially internal, there is no 
need to apply any boundary conditions to them. In the 
numerical solution additional cells are attached to but 
outside these regions as shown in Figure 1. The flow 
parameters in the extra cells are set equal to the flow 
parameters in the corresponding cells inside the computa- 
tional domain which are displaced in the Y direction by 
the cascade blade spacing, S. 

3.3. INLET CONDITIONS 

The U component of velocity, the flow angle 6, the 
gas stagnation enthalpy, H, and the entropy, S, are kept 
constant at the inlet. This leads to a unique definition 
of the flow parameters at the inlet boundary as follows 
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U^ n = const; V^ n /U^ n = tan0 = const ; 



Y P< 
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The explicit formula for V. , p. , and P. are therefore 
* in' in' in 
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3.4. OUTLET CONDITIONS 



At the outlet boundary we define only pressure p out 
and for all other flow parameters apply continuation con- 
ditions. This means that U . , V . , and p^ . are set 

out out *out 

equal to the values of U, V and p one point ahead of the 
outlet boundary. 
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IV. THE NUMERICAL METHOD 



In the program there is a choice of using one of the 
following numerical methods: 

a) First order accurate Godunov Method (Ref. 1) 

b) Higher order extension of the Godunov method (Ref. 2) 
The higher order method is an extension of the Godunov method 
and was used as it is described in Ref. 2. The piecewise- 
linear algorithm was used rather than piecewise-parabolic . 

Here we will describe the Godunov method as it is imple- 
mented by the program. 

Let us assume that the grid is covering the computational 
domain as is shown in Fig. 2(a). We will also consider a 
fixed grid with the indexing of lines as shown in Fig. 2(a). 
The flow parameters related to the center of the cell will 
have a fractional index: j + 1/2, k + 1/2. The cell boun- 
daries will have one fractional and one integer index: 
j, k + 1/2 or m + 1/2, n. The parameters at the time t 
will have subscript index (i.e., Pj + \/2 k + 1/2 ^ < ^ s_ 

tinction from the parameters at the time t + At which will 
have superscript index: + ^ + 

In this notation Eq. (1) is approximated to the first 
order of accuracy by the following system of equations in 
finite differences: 
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a j + ?2 ' k+J2 = a. 



At 
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Here, and in the following text with capital letters 
R, U, V, P, etc., we will note the parameters calculated at 
the edges of the zones. To obtain the value of these param- 
eters we calculate the Riemann problem with left and right 
states in the centers of adjacent cells (Ref. 3). For exam- 
ple, to calculate the flow parameters on the edge j,k+^ 
we assume that the left state is located in j _3 5 ,k+^ and the 
right state in j+h,k+% (see Fig. 2) . The velocity vectors 
(U,V) j_j, and ( U 'V)j+}j fc+ij should be expressed in terms 

of normal and tangent components of velocity to the edge 
(j,k+l) - (j,k). Then the Riemann problem is solved for 
following parameters: 



Left state: (p, P, U. 



n' V t^j-^,k+5s 



( 8 ) 



Right state: (p, P, 



V t* j+h,k+h 



where 



U n - the velocity component normal to the edge 



( j , k+1) - ( j ,k) 



V. - the velocity component tangent to the edge 
( j , k+1) - ( j ,k) 



Solution of this Riemann problem will give: 




Then we transform the component of the velocity back to the 
cartesian coordinate system and will have the final value of 
the parameters on the cell edge: 



(R, P, U, B) 



j /k+*j 



That is a brief description of the Godunov method. 

A more detailed description of the method can be found in 
Reference 4. 

The program includes higher order extension of the Godunov 
method proposed by Collela (Ref. 6) which, in order to im- 
prove accuracy, adds the following steps before calculating 
the Riemann problem: 

1) Compute linear profiles of the dependent variables 
in each zone by interpolating slopes at the centers 
of the zones, subject to certain monotonicity 
constraints. This gives rise to a global distri- 
bution of the dependent variables which is piece- 
wise linear, linear in each zone, with jump dis- 
continuities at the edges of the zones. 

2) Correct left and right states by tracing approxi- 
mate characteristics and solving difference ap- 
proximations to the characteristic equations. 

Only then is the Riemann problem solved for the 
corrected right and left states and the values 
obtained by this solution are assigned to the 
corresponding edge of the grid cell. 
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V. COMPUTER PROGRAM 



The program described here was developed on the NPS 
IBM 370-3033 computer. 

Since the computer code was first developed in Pascal, 
which implies structured programming, the program structure 
is simple. The Data Flow Diagram in Fig. 3 illustrates the 
general flow of data in the program and the subroutine 
functions. The diagram is self-explanatory. 

5.1. GENERAL DESCRIPTION 

The computer program consists of 18 routines. Here we 
will describe, in short, the general function of each routine. 
To understand fully the place of each routine and the program 
structure, the reader should follow the Data Flow Diagram 
in Fig. 3 while reading the written description. 

The subroutine MAIN calls all other subroutines in the pro- 
gram. It controls the execution of the program by checking 
two constants ITRNUM and ITRCH . ITRNUM is the number of 
iterations performed. ITRCH controls the number of iterations 
from the last call to the write routines (WRVELC, WRPRES, 

WRDENS or WRINTG) . Three important constants are defined 
in MAIN: 

MAX = number of grid points in x direction 

MAY = number of grid points in y direction 

LEV = intergration level. If LEV=10 the integration 

will be on two grid levels (not recommended) . If 
LEV=9 it is only one level of integration. 

Subroutine RDDATA reads the grid definition data which is 

created on File 3 by the successive execution of the programs 
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GRID and TRNEW . The program GRID creates a two dimensional 
grid of the H-type for the given cascade definition. The 
program TRNEW redefines the grid to the form needed for 
application of the method and calculates the cosine and sine 
of angles defining the direction normal to the cell edge. 

Program TRNEW also defines the array REG which simulates the 
computational domain by assigning to each physical point in 

the domain an integer number. The description of the programs 
GRID and TRNEW will not be given in the present manual. Here 
we shall assume that File 3, for a given cascade, has been created. 
This file contains the following data: 

REG - array of integer numbers which simulates the 
structure of the computational domain 

Al - array containing Ay values on the cell edges and 

the area of the cell located in the point (j+^k+Jj) 

Bl - array containing Ax values on the cell edges 

TCOS - cosine of the angle between the normal direction to 
the cell edge and direction of the x-axis 

TSIN - sire of the angle between the normal direction to 
the cell edge and direction of the x-axis 

VF - x coordinates of the points on the grid covering 
the computational domain 

VY - y coordinates of the points on the grid 

If the program continues integration from a certain step 
for which the values of the flow parameters are known, these 
values should be stored on File 4 (program creates new data on 
File 4 after a certain number of iterations) . In this case 
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subroutine RDDATA reads three additional arrays from File 4: 

RO - density 

P - pressure 

UB - velocity component in x direction 

VB - velocity component in y direction 

Subroutine RDDATA defines the following constants: 

ITNO - number of iteration in the solution of the 
Riemann problem 

RZ - characteristic density 

PZ - characteristic pressure 

XZ - characteristic length 

G - gamma 

Subroutine DEFINE defines the following parameters : 

ITRNUM - number of iteration 
TIMTOT - total time 

RADL - array, radii of curvature on the lower side of 
the profile 

RADUP - array, radii of curvature on the upper side of 
the profile 

If the program is executing from the beginning, the intial values 
of the array RO, P, UB, VB are defined in the subroutine DEFINE. 
The initial enthalpy - ENTALP, entropy - ENTROP and tanQ - TAU 

are calculated in this subroutine. 

Subroutine TIME calculates the time step for the next 

integration based on CFL criterion -TJ. This subroutine writes 
the values of ITRNUM, TJ, TIMTOT, MOMENT, ans MASS where 
MOMENT - integral of momentum in x direction 
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MASS - integral of mass 

Subroutine INOUT defines the inlet and outlet boundary 
conditions for the flow parameters. In this routine we cal- 

t 

culate the values of the flow parameter in the points (2,k) 
and (MAX-1, k) , where k has the values corresponding only to 
the centers of the cells. 

Subroutine BOND defines boundary conditions for the com- 
putational domain, except inlet and outlet boundary conditions 
which are defined in the routine INOUT. Two kinds of boun- 
dary conditions are calculated in this routine, as a function 
of the cell type indicator REG(j,k): 

1) Boundary condition at the upper and lower surfaces 
of the cascade profile 

2) Symmetry boundary conditions for the regions 1-3, 

5-7, 2-4, and 6-8. To calculate the values of the 
flow parameters at the cell edges this subroutine 
calls routine ITER to solve the Riemann problem at 
the edges. 

Subroutine RIMSOL defines and solves the Riemann problem 
along k lines between the centers of the cells. This routine 
has an option to calculate the more exact values for the left 
and right states calling routine CHARY. A condition state- 
ment in this routine controls call of the routine CHARY that 
allows the use of the same code for the Godunov method and its 
higher order accurate extension. 

Subroutine RIMSLX defines and solves the Riemann problem 
at the edges of the cells in the y line directions. In the 
case of higher order accurate integration, this routine should 
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call the subroutine CHAR. That call is controlled by a con- 
dition statement. 

Subroutine INTGR performs integration of Eq. (.1) in accord- 
ance with the numerical formula (7) . This subroutine calculates 
the values of MOMENT, MASS and ENERGY. 

Subroutine WRVELC calculates and writes the values of 
the local Mach number for the center of each cell on File 6. 

Subroutine WRPRES writes the values of the normalized 
pressure P(j,k) in the centers of the cells on File 6. As 
an option this subroutine calculates and writes the pressure 
coefficient C instead of P(j,k). 

■C'* 

Subroutine WRDENS is included as optional and it writes 
the values of the array RO(j,k) (normalized density) on 
File 6. 

Subroutine UUE writes the values of the local Mach num- 
ber corresponding to the profile surface on File 6. 

Subroutine WRINTG writes complete content of the arrays? 

RO, P, UB AND VB, on File 4. That allows in the next execu- 
tion of the program, beginning integration from the last it- 
eration of the previous execution of the program. 

Subroutine VISCO is optional and is not used in regular 
computation. This routine adds the numerical viscosity of 
Lapidus type (Ref. 5) to the conservation equations. 

Subroutine RIEMAN calculates the Riemann problem for the 
defined left and right states. 
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Subroutine ITER (called only from subroutine RIEMAN) 
calculates pressure on the interface for the solution of the 
Riemann problem. 

Subroutine CHAR corrects values of the left and right 
states solving the equations of gas dynamics in the charac- 
teristics form. In CHAR this correction is done in the 
x direction. 

Subroutine CHARY is the same as CHAR, only for y direction. 

Subroutine MONO, is called only form CHAR or CHARY. This 
subroutine computes linear profiles of the dependent variables 
in each zone by interpolating slopes at the centers of the 
zones, subject to certain monotonicity constraints. 
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VI . INPUT 



All external input to the program is written on File 3. 
Data on File 3 is written in Format (2613) for array REG ( j , k) 
and in Format (5F15.12) for all other arrays. If the user 
wants to continue iteration from the last iteration of the 
previous execution of the program, then the user should input 
the values of the arrays of the flow parameter. This data 
is recorded on File 4 with Format (5F15.12). 

All external input arrays are read from the subroutine 
RDDATA. The inlet parameters of the flow and the outlet pres- 
sure are defined in the subroutine DEFINE from line 1560 to 
1710. 
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VII. OUTPUT 



All output of the program are recorded on File 6 and 
File 4. The following subroutines write output on to 
File 6: 

TIME, WRVELC, WRPRES, WRDENS , UUB . 

Only subroutine WRINTG writes on File 4. 
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VIII. RESULTS 



The Fortran version of the computer code was tested for 
a number of problems for which the solution was obtained with 
the Pascal version. The conclusions from this testing are: 

1) Fortran program is - 5Q% faster than Pascal program, 
when executed with the option Optimize (2) on the 
Fortran H compiler with single precision. 

2) Execution with Fortran G compiler is 30% faster than 
Pascal. 

3) Execution with Fortran H compiler with Optimize (1) 
option is 30% faster than Pascal. 

4) Execution with Fortran H compiler with double pre- 
cision with Optimize (2) option is only 25% faster 
than Pascal. 

5) The differences in results are caused by different 
levels of precision involved in the calculations. 
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2 




Fig. 1. 



Computational grid and boundary of the 
computational domain. 



We look for the solution of J:he_system of equations 
represented by equation (1) in the computational 
domain shown in Fig. 1 for t->-+“ with the 
following conditions at the domain boundaries: 

a) Inflow along segment 1-2 

b) Outflow along segment 7-8 

c) Solid wall along segments 3-5 and 4-6 

d) Cascade : Periodicity between segments 

1-3 and 2-4 and between segments 5-7 and 
6-8 

Channel : Solid wall along segments 1-7 

and 2-8 
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f 



A 

Y 






a) Grid Location in the Coordinate System 





b) Single cell 

Fig. 2. Computational Space Notation. 
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INTGR 




Fig. 3. (continued) 
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